#-------------------------------------------------------------
# GAMLSS(Generalized additive model for lambda,scale and shape)
#-------------------------------------------------------------
library(gamlss)
library(data.table)
library(ggplot2)
library(ggsci)

setwd("D:/projects/Vision")

df <- fread("./Data/prev.csv")
df <- df[!is.na(odstdse),]
df[,newstdse:=odstdse+12.5]
df[,min(newstdse)]
df[,max(newstdse)]
df2 <- df[,.SD[1],by=.(card_no)]
rm(df)

gam <- gamlss(stdse~cs(semster,3), 
              data=df, 
              family=NO
              ) 
save(gam,file="./Data/gam.Rdata")
rm(gam)
# BCCG
gamlms2 <- gamlss(newstdse~cs(semster,3), 
                 sigma.formula = ~semster,
                 methods=CG(),
                 data=df2, family=BCCG)
save(gamlms2,file="./Data/gamlms2.Rdata")

load("./Data/全人群分位数.Rdata")
plot(gamlms2)
wp(gamlms2)
summary(gamlms2)

prediction=centiles.pred(gamlms2,
              type = "centiles",
              xname="semster",
              cent = c(1,5,10,25,50,75,90,95,99),
              calibration = F,
              xvalues = seq(1,18,1)
              )
prediction[,2:ncol(prediction)]=prediction[,2:ncol(prediction)]-12.5
prediction=tidyr::pivot_longer(prediction,2:ncol(prediction))
prediction$name <- as.numeric(prediction$name)

# 全人群
ggplot(aes(x,value,group=name,color=factor(name)),data=prediction)+
  geom_point()+
  geom_smooth()+
  scale_color_manual(name="Percentile",
                     values = c("#D75C4E","#F39764","#C89836","#C8C0DE",
                                "#B5CE4F","#97D2C5","#7EABD0","#52A9DC",
                                "#2D8877")
                     )+
  ylab("SER")+
  scale_x_continuous(name = "Semster",breaks = seq(1,18,1))+
  theme_classic()+
  theme(text = element_text(family = "serif",size=12))
ggsave("./Plots/全人群分位数.jpg",dpi = 300)
save(prediction,file="./Data/全人群分位数.Rdata")
rm(gamlms2,prediction)
rm(df2)
gc()


# 非近视人群
df2 <- df[myopia==0 & odstdse>=-0.5,.SD[1],by=.(card_no)]
gamlms3 <- gamlss(newstdse~cs(semster,3), 
                  sigma.formula = ~semster,
                  methods=CG(),
                  data=df2, family=BCCG)
load("./Data/非近视人群分位数.Rdata")

prediction=centiles.pred(gamlms3,
                         type = "centiles",
                         xname="semster",
                         cent = c(1,5,10,25,50,75,90,95,99),
                         calibration = F,
                         xvalues = seq(1,18,1)
)
prediction[,2:ncol(prediction)]=prediction[,2:ncol(prediction)]-12.5
prediction=tidyr::pivot_longer(prediction,2:ncol(prediction))
prediction$name <- as.numeric(prediction$name)
ggplot(aes(x,value,group=name,color=factor(name)),data=prediction)+
  geom_point()+
  geom_smooth()+
  scale_color_manual(name="Percentile",
                     values = c("#D75C4E","#F39764","#C89836","#C8C0DE",
                                "#B5CE4F","#97D2C5","#7EABD0","#52A9DC",
                                "#2D8877")
  )+
  ylab("SER")+
  scale_x_continuous(name = "Semster",breaks = seq(1,18,1))+
  theme_classic()+
  theme(text = element_text(family = "serif",size=12))
ggsave("./Plots/非近视人群.jpg",dpi=300)

# 近视人群
df2 <- df[myopia==1 & odstdse<(-0.5),.SD[1],by=.(card_no)]
gamlms4 <- gamlss(newstdse~cs(semster,3), 
                  sigma.formula = ~semster,
                  methods=CG(),
                  start.from = gamlms4,
                  control = gamlss.control(n.cyc = 40),
                  data=df2, family=BCCG)
save(gamlms4,file = "./Data/近视人群屈光度.Rdata")
